     function [ymat,inov]=varp(in,arm,insig); 
%   function [ymat,eps]=varp(in,arm,insig); 
%   Simulates an n-dimensional VAR of order P 
%   Inputs
%   -------
%   in.n       dimension 
%   in.lags    Var(P)
%   in.drop    number of observations to drop 
%   in.take    number of observations to take 
%   arm        bphi = [ phi(1) phi(2)    phi(p)] 
%   insig      [n x n]covariance matrix of the errors 
%   Output
%   ------
%   ymat      [in.take x n] matrix of generated series 
%   inov      [in.take x n] matrix of innovations 
%   ---------------------------------------------------
%   Alejandro Justiniano  jalejand@princeton.edu 
%   2/26/03 
%   ---------
	ss=in.lags*in.n; 
	[nch,ssch]=size(arm); 
	% Write in companion form  
	ff=zeros(ss,ss); 
	ff(1:in.n,:)=arm; 
	ff(in.n+1:end,1:ss-in.n)=eye(ss-in.n); 
	% Find eigenvalues to determine dynamics of ths system
    v=abs(eig(ff)); 
	ch1=any( v >= 1); 
	ch2=any(v == 0); 
	chk=ch1 + ch2; 
	if chk > 0; 
        warning('System in non-stationary or contains exact unit roots'); 
	end;
	clear ch1 ch2 matf s v d; 
	if ssch ~= ss; 
        error(' Wrong column dimension for lag matrix '); 
	end; 
	if nch ~= (in.n); 
        error(' Wrong row dimension for the lag matrix '); 
	end; 
	clear nch ssch; 
	chsig=chol(insig); 
	inerr=(randn(1,in.n)*chsig)';
	to=in.take+in.drop; 
	ymat=zeros(to+in.lags,in.n); 
	ymat(in.lags,:)=inerr;
	inov=randn(to+in.lags,in.n)*chsig; 
	% First observation corresponds to a draw from the residuals 
	ii=in.lags+1; 
	while ii<=to+in.lags;   
        xt=ymat(ii-in.lags:ii-1,:); 
        xt=flipud(xt); xt=reshape(xt',in.n*in.lags,1); 
        yt=arm*xt + (inov(ii,:))'; 
        ymat(ii,:)=yt'; 
        ii=ii+1; 
	end; 
	ymat=ymat(in.lags+in.drop+1:end,:); 
	inov=inov(in.lags+in.drop+1:end,:); 